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ABSTRACT 

We present detailed three-dimensional radiative-hydrodynamical models of the well 
known irradiated exoplanet HD189733b. Our model solves the fully compressible 
Navier-Stokes equations coupled to wavelength-dependent radiative transfer through- 
out the entire planetary envelope. We provide detailed comparisons between the exten- 
sive observations of this system and predictions calculated directly from the numerical 
models. The atmospheric dynamics is characterized by supersonic winds that fairly 
efficiently advect energy from the dayside to the nightside. Super-rotating equatorial 
jets form for a wide range of pressures from 10~ 5 to ~ 10 bars while counter rotating 
jets form at higher latitudes. Calculated transit spectrum agree well with the data 
from the infrared to the UV including the strong Rayleigh scattering seen at short 
wavelength, though we slightly under-predict the observations at wavelengths shorter 
then ~ 0.6^m. Our predicted emission spectrum agrees remarkably well at 5.8 and 
8/nm, but slightly over-predicts the em ission at 3.6 and when compared to the 

latest analysis bv lKnutson et all (|2012f ). Our simulated IRAC phasecurves agree fairly 
well with the amplitudes of variations, shape, and phases of minimum and maximum 
flux. However, we over-predict the peak amplitude at 3.6/ito and 4.5^m, and slightly 
under-predict the location of the phasecurve maximum and minimum. These simula- 
tions include, for the first time in a multi-dimensional simulation, a strong Rayleigh 
scattering component to the absorption opacity, necessary to explain observations in 
the optical and UV. The agreement between our models and observations suggest that 
including the effects of condensates in simulations as the dominant form of opacity 
will be very important in future models. 

Key words: hydrodynamics, radiative transfer, methods: numerical, planets and 
satellites: atmospheres, planets and satellites: HD189733b 



1 INTRODUCTION 

Among the ~ 630 known exoplanets, HD189733b currently 
boasts the most numerous and comprehensive set of obser- 
vations to date. A member of the well know class of short- 
period, highly irradiated gaseous planets (hot Jupiters) its 
relative proximity to earth and the favorable star to planet 
radius ratio have made it the preferred target for many 
different groups studying radial velocity, primary and sec- 
ondary transits, transit spectra, emission spectra, variabil- 
ity, and phasecurves across a wide range of wavelengths. 
Given this wealth of knowledge, it is the best candidate 
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for detailed comparisons between model predictions and ob- 
servations. In this paper we compare the results from our 
three-dimensional radiative-hydrodynamical simulations of 
HD189733b with this wide array of observations. 

Hot- Jupiters have been found around approximately 
10% of solar-type stars. Their short orbital periods (~ 1 — 4 
days) suggests these planets are subject to strong tidal forces 
that synchronize and circularize their orbits. As a result one 
side perpetually faces its host star, receiving intense irra- 
diation, while the other side receives no stellar irradiation. 
In the absence of any hydrodynamics, the dayside reaches 
temperatures of several thousand degrees, while the night- 
side temperature, maintained by the internal energy of the 
planet, would remain at several hundred degrees. However, 
this large temperature differential across the planet acts as 
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an enormous driving force for the atmospheric gas generat- 
ing supersonic winds that advect energy across the planet. 
The resulting temperature distribution depends crucially 
on the non-linear coupling between the radiative transfer 
of energy and hydrodynamic transfer of energy. The spa- 
tially varying temperature across the planet leads to vary- 
ing chemical composition, and emission and absorption ef- 
ficiencies and is thus directly tied to observable properties. 
HD189733b is a fairly typical example of such a planet. 



There are a wide range of numerical approaches 
that have been taken to address the coupled radia- 
tion hydrodynamics of irradiated gas giants. Though cou- 
pled, it is convenient to separate a given studies ap- 
proach for the radiative transfer from that for the hy- 
drodynamics. Though not an exhaustive list, previous ap- 
proaches to radiative transfer include relaxation meth- 



ods (i.e. Newtonian heating) (IShowman & Guillot 


2002 


Showman et al. 20081; ICoooer & Showmanl 120051. 


2ooe 


Laneton & Laughlinl 2007. 20081: Menou & Rauscherl 


2009; 


Rauscher & Menou 2010: Pcrna et al.ll201dl). kinematic con- 



straints de s igned to 
120031 . 120081 : IRauscher et al 
sion (FLD) with (|ppbbs-Dixon et al 
|Burkert et all |200dF 



120081). 3D flux -limited diffu- 
ixon et al.l |2010|) and without 
Dobbs-Dixon fc Linl 2008) a separate 



radiation component , dual grey ID two-stream appr ox- 
imation jHeng et all l201ll ; IRauscher fc Menoul l2012al lbl). 
3D FLD coupled to waveleng th dependent stellar irradia- 
tion (|Dobbs-Dixon et al.ll2012p and ID frequency -dependent 
radial radiative transfer ( Showman et al.l [2009). Previous 
methods used to solve the hydrodynamical portion in- 
clude s olving t he equivalent bar o tropic e quations llCho et al. 
2001 I 2008T IRauscher et~ai1 l200Sl : IRauscher fc Menoul 
2012al lbJ), th e shallow water equations (Lang ton fc La ughlin 
20071 , l2008h. the primi t ive eq u ation s IShowma n fc Guillotl 
20021 ; IShowman et al] l200Sl. 120091; ICooper fc Showmanl 
Menou fc Rauscherl I2OO9I; iPerna et al] l20ia 



Rauscher fc Menoul I2OI0I: iHeng et al.l l201ll), Euler's equa- 



tions ijBurkert et aD 120051 : iBobbs-Dixon fc Linll200Sl ) , and 
the Navier-Stokes equations ( Dobbs-Dixon et al] 120101 . 
2012). There are a wide range of assumptions in all these 
approaches. Currently the most sophisticated approach to 
hydrodynamics utilizes the Navier-Stokes equations, while 
that for radiative transfer are ID frequency-dependent mod- 
els. 



2 MODELING METHODOLOGY 

In this paper we present results from a numerical model 
coupling the fully compressible 3D Navier-Stokes equations 
to a two-stream, frequency dependent radiative transfer cal- 
culation. Though we solve the Navier -Stokes equations in 
a manner similar to previous p apers l|Dobbs-Dixon fc Linl 
120081 ; iDobbs-Dixon et al] |2010L |2012| ) the radiative trans- 
fer methodology is presented here for the first time. Ob- 
servables (phasecurves, spectra, etc.) are calculated utilizing 
the results of the coupled radiative-hydrodynamical models. 
We describe all three components (hydrodynamics, radiative 
transfer, and the calculation of observables) in detail below. 



2.1 Hydrodynamics 

The hydrodynamical portion of the model solves the 
fully compressible three-dimensional Navier-Stokes equa- 
tions given by 
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representing the momentum, continuity, and thermal en- 
ergy equations respectively. The three-dimensional velocity 
is given by u, p and P are the gas density and pressure, and 
Q is the rotation frequency. The gravitational acceleration, 
g is taken to be constant and purely radial, v = 77/p is the 
constant kinematic viscosity. We have neglected the coef- 
ficient for the bulk viscosity. The thermal energy equation 
contains heating and cooling contributions from the radial 
gradient of the radiative flux Fr, the incident stellar energy 

and the viscous heating D v . The formalism for Fr and 
S+ is detailed in Section (|2.2[) and that for D v can be found 
in lKlev fc Henslerl l| 19871) . 

We solve these equations in spherical coordinates with 
a resolution of (N r ,N^,Ne) = (100, 160,64) where <f> is the 
longitude and 9 the latitude. At the equator this corresponds 
to cells that are 2.25° by 3°. Flo w over the pole is ac c ounte d 
for via the method described in lDobbs-Dixon et al] |2012). 



Given the wide range of approaches utilizing a wide 
range of numerical codes, HD189733b provides an excellent 
target for not only testing how well these models perform but 
also exploring the physical processes included. Therefore, in 
this paper we present three-dimensional models specifically 
tuned for HD 189733b utilizing a model coupling wavelength 
dependent two-stream approximation, wavelength depen- 
dent stellar irradiation, and the fully compressible Navier- 
Stokes equations. Section (2) presents the details on our 
modeling methodology including both the dynamics and ra- 
diation in the simulations and calculating observable signa- 
tures utilizing the models results. Section (3) presents results 
illustrating both the temperature and dynamical structure 
across the planet and detailed comparisons to observations. 
Section (4) concludes with a discussion of our results. 



2.2 Radiative Transfer 

In highly irradiated planets, transfer of energy via radiation 
plays a very important role. For HD 189733b this radiation 
can be broken into two distinct contributions: the irradiation 
the planet receives directly its host star and the local radia- 
tion comprised of re-radiated stellar energy and the flux from 
the planets interior. Though we include a full wavelength de- 
pendent treatment, the incoming stellar irradiation can be 
largely characterized as short wavelength (visible) and the 
re-radiated characterized as long wavelength (infrared). 

2.2.1 Two-Stream Approximation 

To address radiative transfer of the re-radiated energy we 
have developed a frequency dependent, two-stream approxi- 
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mat ion (|Mihalaslll97sh for the radial radiative flux, separat- 
ing it into multiple upward (i*t,fc) an d downward (i*4.,b) prop- 
agating channels. The governing equations in each wave- 
length bin can be written as 



dFA 



dF ub 
dn 



= -F Ub + S b , 



(4) 
(5) 



We have subdivided the full spectrum into 3 wave - 
length bins (identi c al to those in IShowman et al.l (|2009h ; 
iDobbs-Dixon et all HoH)). Assuming the gas emits as a 
blackbody, the source function Sb in Equations Q and Q 
in a given bin can be written as 

pv 2 [b] 

S b = 7T / B v (T, v) dv (6) 

J V1 [b] 

These equations require two boundary conditions. For 
these we choose Fj.^ (r = R p ) = at the surface and 
Ff,b { r = Rbot) = Sb (Tbase) at the interior, where we fix the 
total upward flux 



y ] Sb (rbot 



dT 
dr 



Ar Ar = F int 



(7) 



to the value required to match the observed radius. Note 
that though we set the downward flux to be zero at r = R p 
in Equation ([5]), stellar heating is accounted for directly in 
the thermal energy equation (Equation p|). This term is 
discussed further in the next subsection. The optical depth 
in Equations Q and ([5} are calculated utilizing the averaged 
opacity within that bin and a diffusivity factor, 



Tb (r) = —1.66 / pK b dl 



(8) 



The diffusivity factor of 1.66 is an approximation that ac- 
counts for an exponential integral that arises when tak- 
ing the first m oment of the intensity to calculate the flux 
(Elsasser 1942). The opacity Kb is discussed in more detail 
below. 

We solve Equations Q and (J5]) in the standard way, uti- 
lizing an integrating factor e T . Together with our boundary 
conditions we find 



Fl,b (Tb) 



S b e 



(9) 



Note that if we wished to include the stellar heating directly 
into these equations there would be an additional exponen- 
tially attenuated term Sb,*e~ Tb added to the above equation 
for with 



(— 7T / B u (T*,v)dv 



(10) 



The equation for the upward flux does include an initial 
source term due to the flux from the interior, and can be 
written as 

F t ,i {n) = s b (T bot ) e -K^-0+ j b s b e-«- Tb )dri(ii) 

•> T b,bot 



Finally, the net radial flux (defined with positive indicating 
outward flux) can be computed as 



(n) = ( F t,b (n) - F Ub (r b )) . 

b 



(12) 



However, it is important to note that Equation (1121) 
only works in the upper at mosphere. As demonstrated by 
iRauscher fc Menoul (|2012ah . in regions below the photo- 
sphere the temperature-pressure profile becomes isother- 
mal. As the density increases, the integrated optical depth 
across a grid cell becomes very large, and the profile be- 
comes dependent on the numerical resolution. The resolu- 
tion required to compute an accurate temperature profile 
at depth using Equation (|12[) far exceeds computational 
capabilities. Results from other grou ps that utilize some 
form of two-stre am approximation |Showman et al l 120091 ; 
iHeng et alJl201ll ) seem to contain similar isothermal regions 
at depth, though they do not discuss this is s ue. To sur- 
mount this issue, we follow IRauscher fc Menoul (|2012al ) and 
transition to a diffusive scheme at r = 2.5. This is done 
in each wavelength bin independently given that the pho- 
tosphere can be at significantly different depths across the 
spectrum. Given the larger optical depths, a diffusion ap- 
proximation is entirely adequate in the deeper atmosphere. 
Below the active weather layer, radial radiative diffusion re- 
mains the dominant form of energy transport all the way 
do wn to the radiative- c onvec tive boundary. Finally, as seen 
in IDobbs-Dixon et al.l (|2010T ). non-radial radiative transfer 
does play a role in influencing the energetics. Though incon- 
sequential near the sub-stellar point, the horizontal flux of 
radiation can play an important role. We have neglected the 
non-radial contribution here. 



2.2.2 Stellar Heating 

The irradiation of the planets atmosphere from the central 
star is directly accounted for in Equation Q for the ther- 
mal energy of the gas. We choose this procedure rather then 
inserting an additional boundary term into Equation (O for 
Fi t b as it more accurately captures the slant-optical depth 
associated with the path of stellar photons through the at- 
mosph ere. The formalism here is identical to the procedure 
used in IDobbs-Dixon et alj l|2012T ). The stellar heating term 
in the thermal energy equation can be written as 



S* = 



2 nb 
- — - dr b 



dr 



nB v (T*,v)dv, (13) 



"6,1 



where B v (T*, v) is the stellar blackbody for HD 189733 A 
taken from the Kurucz modelfl _R* and a represent the 
stellar radius and semi-major axis. The slant-optical depth 
is accounted for by including the cosine of the angle 
between the normal and the incident stellar photons. The 
optical depth to stellar photons is calculated as 



Tf,.., 



pKb^dl. 



(14) 



2.2.3 Opacities 

The details of the atmospheric opacity is a crucial param- 
eter in determining the energy deposition and re-radiation, 
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subsequently playi ng an important role in the overall atmo- 
spheric dynamics. iDobbs-Dixon fc Linl (|2008T ) explored the 
role of changing the opacity within the context of a grey 
flux-limited diffusion (FLD) approximation. The results of 
this study indicated that as opacity of the atmosphere is 
decreased, stellar energy is deposited at larger pressures, 
where the dynamics is more effective at redistributing it to 
the nightside. In this paper we are utilizing a much more so- 
phisticated radiative transfer routine, but the general prin- 
ciple still holds, increasing opacities causes the energy to be 
deposited much higher in the atmosphere. 

Here we are utilizin g the f requency dependent opacities 
k„ of Sh arp fc Burrows] (|2007l ) which assume solar compo- 
sition gas. The dominant molecules relevant for HD 189733b 
in these calculations are water, methane and carbon monox- 
ide. The frequency dependent opacities are averaged within 
each wavelength bin utilizing a Planck mean, 



Kv {p,T)B v (T)du 



(15) 



However, we find that the exact form of the averaging 
method does not have a significant effect on the radiative 
solution. For the re-radiated component of radiation (Equa- 
tions ([9]) and (|lip ) the local gas temperature is used in 
B v (T). When computing the contribution from the incident 
stellar irradiation (Equation (|14|) ) we replace this with the 
stellar blackbody, B u (T*). 

After significant experimentation, we found that 
radiative-hydrodyn amical models util i zing s imply the opac- 
ities directly from ISharp fc Burrows] |2007) did not agree 
well wit h the observations . Comparison to emission 
spectra jDeming et all 120061; IKnutson et al] |2007t 120091 ; 



ICharbonneau et al* 20081 ; Knutson et al. 20121 ) indicated 
these initial models were too cool in t he upper atmo- 
sphere, and comparison to transit spectra (|Pont et al.ll2008l ; 
ISing et alj|201ll ; iPont et alj|2012l) indicated that the exist- 
ing Rayleigh contribution in the opacities was woefully in- 
adequate to explain the strongly varying transit radius with 
decreasing wavelength. To address this issue we have sup- 
plemented the opacities with two addition components, a 
uniform grey component K grey at all wavelengths and an 
component that scales as A -4 , as would be expected from 
Rayleigh scattering. This extra opacity can be written as 



(A) — Kgrey + ^RS 



A 



0.9/im 



(16) 



We have run a number of simulations varying the magni- 
tude of these two components. Setting K gr ey = 0.035cm 2 /g 
and krs = 0.6cm 2 j g provides the best fits to the observa- 
tions. Note that the magnitude of the grey component is 
th e same as that requi red in the one-dimensional models 
m IKnutson etHI l|2012h . This extra opacity is added self- 
consistently to all the opacities utilized in the code and in 
the post-processing described in the next sub-section. This is 
the first multi-dimensional radiative hydrodynamical model 
to include such a strong Rayleigh scattering component. 
Note however that we are treating k b solely as an absorp- 
tive opacity, ignoring its scattering properties. This should 
be included in future work. 



2.3 Calculating Observables from 3D Models 

To determine the transmission spectrum from the hydro- 
dynamical models, we calculate the wavelength dependent 
absorption of stellar light traversing through the limb of the 
planet. This allows us to determine the effective radius of 
the planet and a fractional reduction of stellar flux F*. Ne- 
glecting limb-darkening of the star, this can be expressed 
as 

r(b, a ,\)\ bdMa 

-nt 2 ■ ^ 



intransit 



1(1 



where r (b, a, A) is the total optical depth along a given chord 
with impact parameter b and polar angle a, defined on the 
observed planetary disk during transit. 

To calculate the emission spectrum and phasecurves, we 
integrate inwards along rays parallel to the line of sight. The 
emerging Blackbody flux from the planet at each point on 
the surface is given by 

2tt/ic 2 /A 5 



B(x,y,X) 



~ T dT. 



(18) 



To calculate the apparent dayside luminosity, we integrate 
over the observed disk 



L P (A) = / I P dA, 



(19) 



where dA is taken over the observer's plane and the emerging 
intensity is given by 



A 

Ip(x,y,\) = — I B x (x,y,T)e~ rx dT X . 



(20) 



For both transmission and emission calculations, the density 
and temperature needed to calculate r at each location are 
interpolated from the values in the 3D models. 



3 RESULTS 

In this section we present the results of our 3D radia- 
tive hydrodynamical simulations of the irradiated planet 
HD189733b. We first discuss the wind and temperature 
structures throughout the atmosphere and then compare 
synthetic observations calculated using these results to a 
wide variety of observations. 



3.1 Temperature and Wind Structure 

As discussed in Section (1) the synchronous rotation and in- 
tense irradiation of the dayside by the host star HD 189733 A 
drives winds throughout the upper portions of HD189733b. 
These winds advect energy throughout the atmosphere re- 
sulting in significantly different temperature distributions 
then expected from purely radiative calculations. As is to be 
expected, the details of the dynamics depends on the input 
parameters; Table {TJ lists a number of these parameters. 

The radial distribution of the energy deposited in the 
atmosphere from its host star extends from the upper re- 
gions of the atmosphere (10 bars) down to ~ 10 bars due 
to the complicated wavelength dependent opacities. Cou- 
pled with radial mixing, this leads to varying temperature 
distributions at different pressure levels within the atmo- 
sphere. Figure {TJ illustrates the temperature and zonal (i.e. 
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Parameter 


Value 


F*roti Porb 


2.22 days 


R* 


O.76R 


Rp 


1.138-R.jup 


M p 


1.144M Jup 


a 


0.031AU 


V 


10 7 cm 2 s -1 



Table 1. Physical parameters chosen for simulati ons of 
HD18 9733b presented in this paper taken from I Torres et al.l 
j2008tl- 



longitudinal) velocity at constant pressure surfaces ranging 
from 10~ 4 to 1 bars. There are a number of features of in- 
terest in these plots. The hottest point in the atmosphere 
is near the sub-stellar point just above ~ 0.1 bars, corre- 
sponding to the location where the majority of the stellar 
energy is deposited. The coolest points are high in the atmo- 
sphere in the nightside gyres found at mid-latitudes. Zonal 
jets, super-rotating at the equator and counter-rotating at 
mid latitudes are also quite evident. Jet velocities increase 
with decreasing atmospheric pressure, reaching 6 km/s at 
pressures of 10 -4 bars. The advective signature of the jets 
is seen near the terminators (90° and 270°) as tongues of 
hotter regions stretching onto the nightside. Deeper in the 
atmosphere, where the radiative time-scale is longer, advec- 
tion is more efficient, equalizing temperature and moving 
the hottest point significantly downwind to the east of the 
sub-stellar point. 

Figure © illustrates the behavior of the temperature 
and the zonal velocity with pressure at both the equator 
and mid-latitudes. Several separate temperature peaks are 
seen in these plots. The temperature inversion seen just be- 
low ~ 10 -2 bars near the sub-solar point again corresponds 
with the maximum energy deposition associated with the 
stellar heating. The larger temperatures at longitudes be- 
tween 90 and 100° and pressures between 1 and 10 _1 bars 
corresponds to the heating caused when the westward mid- 
latitude jet deviates toward the equator and interacts with 
the eastward equatorial jet. At mid-latitude the heating be- 
tween 1 and 10 _1 bars is generated by the compression and 
shock-heating associated with converging flows. Much of this 
behavior can also be seen in the zonal velocity shown in Fig- 
ure (0. Deeper in the planet temperatures converge onto a 
single curve that will continue to steepen until becoming 
unstable at the radiative-convective boundary. 

The zonal velocity is almost entirely eastward (super- 
rotating) at the equator. The largest velocities are found 
high in the atmosphere on the nightside of the planet just 
past the eastern terminator. These flows have very high 
Mach numbers, reaching up to 3.5 at the western termi- 
nator where the flow has cooled radiatively but still re- 
tains high velocities. With such high velocities it is likely 
that shocks play an important role in limiting the veloc- 
ity. Though our numerical scheme can capture larger scale 
shocks, it is likely there are additional processes occurring 
on size scales smaller then our grid that contribute to the 
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Figure 1. Temperature in Kelvin (left) and zonal speed in km/s 
(right) at 10~ 4 , 10~ 3 , 10~ 2 , 10 _1 , and 1 bars from top bottom. 
The images are projected onto a cos(longitude) map to highlight 
the importance of the equatorial features in determining observ- 
able properties. The sub-stellar point is located at 0° latitude and 
the equator runs horizontally through the center of each plot. 



overall drag of the flow and that these velocities represent 
upper limits to the actual flow speeds. At mid-latitude, flow 
is in both the eastward and westward directions at low pres- 
sure. As the pressure increases, circumplanetary westward 
jets develop. 



3.2 Comparison to Observations 

HD198733b is currently the most observed extrasolar planet 
given its relative proximity and favorable planet-star ratio. 
Therefore in this section we perform detailed comparisons 
between these observations and synthetic observations cal- 
culated utilizing our numerical model. The methods for cal- 
culating observable metrics from the simulations were pre- 
sented in Section (|2.3p . 

In Figure ((3} we present the calculated transit spectrum 
from the UV through the infrared. The calculated transit 
spectrum is dominated by Rayleigh scattering below ~ l[im 
with molecular features becoming important at longer wave- 
lengths. However, the uniform component of n e that we have 
added to the opacity suppresses the strength of many of 
these features. The observational data shown in Figure l|3"|l 
is a re analysis of data from a number of groups bv lPont et al.l 
l|2012h . Absorption depths reported by a number of groups 



20081; ISing et al.ll2009l. 20 111; 


Gibson et al.ll2012l; Desert et al. 


20 111, 120091; iKnutson et al. 


2007. 12009: lEhrenreich et al. 


20071). Desert et al.l (201lh and Pont et all d2012h suggest 
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P (bars) P (bars) 

Figure 2. Temperature versus pressure (left panels) and zonal 
velocity versus pressure (right panels) at the equator and mid- 
latitudes. The inset color wheel indicates the longitude of a given 
profile. The star is located t the right of the inset color wheel, 
thus the black line denotes the sub-stellar longitude. The eastward 
(positive) direction is counter clockwise. 

recently reanalyzed much of this data utilizing a consistent 
treatment for stellar spot corrections, transit properties, and 
stellar-limb darkening. Our model agrees fairly well with 
these re-analyzed data. We do see deviations at very short 
wavelengths (< 0.6pm) in the Rayleigh scattering portion 
of the spectrum. This is likely due to either a slight temper- 
ature inversion in the uppermost portion of the atmosphere 
that our model is not capturing or a change in grain size with 
height, physics not included in our current model. The very 
uppermost radius of our model corresponds to an absorption 
depth of ~ 2.45% and pressures of several microbars. We at- 
tempt to account for any additional absorption by adding an 
isothermal region above this pressure with the temperature 
within each radial column taken from the simulated region. 
The limits on the radial extent is a numerical constraint as- 
sociated with very low density in this region and it may play 
a role in the deviations we see at wavelengths < 0.6/im. 

The next diagnostic we examine is the emerging spec- 
trum on the dayside of the planet. Figure Q illustrates the 
calculated emission spectrum during the secondary eclipse. 
Again, we find general agreement throughout the infrared, 
slightly over -predicting the n ewest measurements at 3.6 and 
4.5pm from iKnutson et al.l ((2012) , but doi ng very well at 
5.8 an d 8pm. We under-predict the flux from lKnutson et al.l 
(2009) MIPS measurement at 24^m. W e do n ot reproduce 
the detailed shape seen by I Swain et al.l (|2009l ) at ~ 2pm. 
As our overall goal is to utilize observations to validate our 
radiative-hydrodynamical model we plot the actual emerg- 
ing intensity as \L\ in Figure From this figure it is 
clear that the vast majority of the radiative energetics is 
occurring in the region from w 1- 10pm. Outside this win- 
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1 10 
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Figure 3. The transit spectrum for the simulat e d pla net 
HD198733b. Over-plotted are data from iPont et al.l ll2012h as 
solid points, and the band averaged mode values as red trian- 
gles. Bands are taken from Table 6 of IPont et al ] l|2012h . Both 
the data and model exhibit strong Rayleigh scattering at short 
wavelengths and muted molecular features at longer wavelengths. 

dow the emerging intensity is orders of magnitude lower. 
The dynamics (jet formation, advection, meridional circula- 
tion, etc.) will overwhelmingly be driven by radiative energy 
transfer in this window. Therefore, the agreement between 
our model and observations in the IRAC bandpasses gives 
us confidence that we are capturing a majority of the radia- 
tive energetics, despite some disagreement at longer wave- 
lengths. To get a better sense of the planetary emission in 
the IRAC bandpasses, we plot the spatially resolved emis- 
sion from the planet in Figure ((6| during the center of the 
secondary eclipse. These figures illustrate a number of fea- 
tures including the offset of the hot-spot, the advection of 
energy by the super-rotating equatorial jet, and the vary- 
ing level of asymmetry in different bandpasses, indicative 
of the varying pressure levels they probe. The IRAC 4 im- 
ag e compares favorabl y to the 2D seconda ry eclipse map 
of iMaieau et all (|2012l ); Ide Wit et all (|2012T l utilizing 8pm 
data. 

One of the most important diagnostics of atmospheric 
dynamics on extrasolar planets are multi-wavelength phase- 
curves of the planet throughout its entire orbit. As vary- 
ing phases of the planet come into view, we are able to 
directly compare the observed flux to that predicted by 
models. Deviations from radiative equilibrium indicates ad- 
vective contributions to energy redistribution via the at- 
mospheric dynamics. The presence of a strong circumplan- 
etary jet wa s first robustly dem onstrated for the planet 
HD209458b (|Knutson et al] l2007h . Therefore, we examine 
the multi-wavelength phasecurves. The amplitude, shape, 
and the phases of the the minimum and maximum of the 
flux are all important diagnostics. In Figure([7| we show 
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Figure 4. The dayside emission spectrum for HD198733b 
shown as the black curve. Over-pl o tted are data from 
Deming etaLI J2006I); ISwain et all ll200Sf); IKnutson et ail j200^ ; 
Charbonneau et alj (|2008h : iKnutson et al.l [ |2009l . l2012n 
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Figure 5. Apparent dayside luminosity AL^ during secondary 
transit clearly illustrating the wavelengths where a vast majority 
of the radiative energy that ultimately drives the dynamics is 
found. 
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Figure 6. The resolved emerging intensity given by Equation i20l 
with units of ergscm~ 1 str~ 1 cm~ 2 s~ 1 , during secondary eclipse 
as would it would be seen in the four IRAC bands across the 
entire face of the planet. The sub-stellar point lies in the center 
of each image. The importance of the circumplanetary jet and the 
varying degrees of asymmetry are evident. 



the simulated phasecurves in the I RAC and MIPS band - 
passes compared t o observations fr om Knutson et al.l (2009); 
lAgol et all (|2009l ); IKnutson et all i|2012l ). In general we find 
remarkable agreement between our models an d the data. 
We do very well in matching the observations o f lAeol et al.l 
(2009) at 8/im, including its distinct, non-symmetric shape. 
At 3.6/j.m we slightly over-predict the flux at secondary 
eclipse (orbital phase of 0.5), and very slightly under-predict 
the flux from the nightside. At 4.5/xm we also over-predict 
the flux at secondary eclipse, but match the nightside flux 
fairly well. Finally we significantly under-predict the flux at 
24/^m. We saw this behavior in Figure (@J) as well, but as dis- 
cussed above the results at 24/im likely have little influence 
on the overall atmospheric dynamics. However, in compar- 
ing our mode l to the pha s ecurve observations it is important 
to not e that I Pont et al.l (|2012h claims that IKnutson et al.l 
(2012) overstates the ability of Gaussian Processing to cor- 
rect for variability on the 1-2 day time-scale. Thus the errors 
in their phasecurve amplitude are likely ~ 5 times too small. 
This suggests our amplitudes may consistent with the data 
given these larger uncertainties. Also seen in Figure© is 
the distinct offset of the phase of maximum and minimum 
fluxes from 0.5 and 0, respectively. Our models agree very 
well with the shape of the curve in the IRAC wavelengths of 
3.6/im, 4.5^m and 8fim. The observations at 24/im are too 
sparsely sampled to definitively compare the phase-offsets. 
Table (|3.2|) lists the observed and simulated phase differ- 
ence from transit and secondary eclipse for the minimum 
and maximum fluxes respectively. Comparison between the 
observations and the models shows that we systematically 
under-predict the magnitude of the offset at all wavelengths, 
save an over-prediction of the minimum at 4.5pm. Taken at 
face-value this suggests that the jet strength may be some- 
what underestimated in our models. 
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Table 2. Comparison between the calculated and observed offsets of minimum and maximum 
fluxes (in hours) from the inferior and superi or conjunction, r espec tively, in the 1RAC and MIPS 
bandpasses. Observed values are taken from [Knutson et al.l ||2012| ) Table 2. Note that duration 
and accuracy of the Sfim and 24/im data was not sufficient define the minimum of the phasecurve 
in these wavebands. 
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Figure 8. The effective time-offset as a function of wavelength. 
Timing variations are primarily due to temperature differences 
between the eastern and western hemispheres. 



Figure 7. The simulated IRAC phasecurves f or HP 189733b 
Over- plotted ar e observed phas ecurves from 
2009) at 24/tm, lAgol et al.l ||2009|) at 8^m, and 
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20121 1 at 4.5 and 3.6A*m. 



Though pha securves provide an excell ent diagnostic of 
the jet dvnamics. lDobbs-Dixon et al.l (|2012T ) suggested a new 
observational diagnostic to determine the strength of the 
jet. In the presence of a strong jet the temperature at the 
eastern terminator (in the direction of the jet) will be some- 
what larger then that at the western terminator. The differ- 
ential scale-heights between the two terminators will alter 
the timing of ingress and egress at a measurable level. This 
effect will change strength at different wavelengths as the 
day/night temperature differential depends on pressure and 
may be detectable with simultaneous observations at multi- 
ple wavelengths. In Figure (JS]) we present the apparent offset 
in the time of central transit when the lightcurves computed 
from our mo del are fit with a sph erical planetary transit 
model (from iMandel fc Arol (|2002h l. We predict that the 
magnitude of this effect will be most visible in observations 
comparing the 3.6 and 4.5/im regions of the spectra (with 



timing differences reaching up to 3 seconds) , and also possi- 
bly at short wavelengths. 

The final diagnostic we explore is the stability of the at- 
mospheric dynamics and thu s the tempe r ature distribution 
throughout the atmosphere. lAgol etall l|2010l ) monitored 
secondary eclipse depths of 7 events spread over 270 days. 
They find very little change in the observed depths, putting 
an upper limit on the variability of the dayside flux of 2.7%. 
To test this, we have calculated the secondary eclipse depth 
variation over 30 orbits. Note that we don't actually have 
to simulate 30 orbits to derive this quantity. Because we 
don't have observational constraints for our models we have 
simply continuously monitored the expected depth over this 
period. We find that the dayside emission from our models 
is extremely stable, with RMS deviations of 0.1%, 0.07%, 
0.07%, 0.07%, and 0.06% at 3.6/oti 4.5/im, 5.8^m, 8>m, 
and 24/im respectively. 

In summary, our model shows fairly good agreement 
with all the observations. It is possible that the jet is some- 
what stronger then we predict and thus more efficient at 
cooling the dayside and heating the nightside, as indicated 
primarily by the 3.6^m phasecurve. In addition, some frac- 
tion of the energy we see emitted at shorter wavelengths 
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Figure 9. The simulated phase curves for the four IRAC bands 
and the MIPS band for simulations with a varying strength of the 
extra grey absorber re e - 



is actually escaping at longer wavelengths. Given the re- 
sults presented in Figure ((5| only a very small fraction 
of the energy at short wavelengths must be transferred to 
longer wavelengths to explain the observations. There are 
a number of possible reasons that we are finding a discrep- 
ancy between the model and the observations. One likely 
culprit is the composition and thus the opacity of the at- 
mospheric gas. As mentioned we are utilizing sola r com - 
position, molecular opacities from lSharp fc Burrows! (|2007T ) 
augmented with a extra wavelength dependent opacity of 

k c = 0.035 + 0.6 ( 015^7) ~ 4 - In Figure© we illustrate the 
effect of varying the strength of a simple grey component 
(Ke) on the phasecurves. Indeed it is clear that increasing 
K e causes the overall emission to shift from the 3.6/im band 
to longer wavelengths. A small decrease at 3.6/xm is more 
then sufficient to cause large shifts at other wavelengths. 



4 DISCUSSION 

We have presented a new set of radiative hydrodynamical 
simulations for the irradiated gas giant planet HD189733b. 
Simulations utilize a frequency dependent two-stream ap- 
proximation to calculate radiative transfer coupled to a 3D, 
fully compressible solution of the Navier-Stokes equations. 
Opacities, crucial for understanding the deposition of stel- 
lar energy and the subsequent re-radiation and cooling of 
the atmosphere, are calculated assuming contributions from 
solar-composition molecules, a wavelength independent grey 
component, and strong Rayleigh scattering. The latter two 
are assumed to come from the presence of grains in the at- 
mosphere. 

The atmospheric dynamics is characterized by super- 
sonic winds that fairly efficiently advect energy from the 



dayside to the nightside. Super-rotating equatorial jets form 
for a wide range of pressures from 10 -5 bars down to ~ 10 
bars. Counter rotating jets form at higher latitudes flanking 
the equator setting up a pattern of three jets from north 
to south pole across the nightside. The westward, counter 
rotating jets become slower at depth, but are also able to 
circumnavigate the planet below pressure of ~ 10 -2 bars. 
Not only is the temperature structure modified by advec- 
tion of energy in these jets, but the interaction between the 
jets also leads to compressional and shock heating, further 
modifying the temperature distribution across the planet. 

We performed detailed comparisons between our model 
and a number of observations including transit spectra, 
emission spectra, eclipse maps, phasecurves, and variability. 
Model transit spectra agree fairly well with the data from the 
infrared to the UV, though they under-predict the observa- 
tions at wavelengths less then ~ 0.6/im. This may be due to 
an increase in temperature or a change in grain size, neither 
of which is captured by our models. Our predicted emission 
spectrum agrees remarkably well at 5.8/im and 8fim, but 
slightly over-predicts the emission at 3.6/im and 4.5/xm whe n 
compared to the latest analysis by iKnutson et all l|2012Ti . 
Our measurement of the phasecurve at 8fim agrees very 
well, reproducing both the peak amplitude and the phase of 
maximum flux. As with the emission spectrum, we over pre- 
dict the peak amplitude at 3.6/im and 4.5/im, though agree 
fairly well with the phase offsets of minimum and maximum 
fluxes. 

Among the other groups studying multi-dimensional 
radiation h ydrodynamics o f irra di ated hot Jup iter at- 
mospheres^ Sjwwmajy^ff] (120091 ). iFortnev et al.l (|2010T ). 
and IKnutson et al.l ( 2012! ') have also performed detailed 
comparisons to HD18 9733b data utilizing the models of 
Sh owman et al.l (|2009f i . As both their model and opacities 
differ from ours it is not unexpected that we differ in some 
of the details. However, despite these differences we do ap- 
pear to present a broadly consistent picture of the dynamics. 
Perhaps the most important difference between our results, 
with respect to the dynamics, is the location of the hotspot, 
indicative of the strength of the w inds and the efficiency of 
advection. Sho wman et al.l (|2009h significantly over-predict 
this efficiency, finding the location of the hottest point some- 
what further downwind then observed. As a result of in- 
creased advection the nightside temperatures they predict 
are larger then observations. To address the discrepancy be- 
tween the observed ph a se of the hotspot and that of the 
model, IShowman et all (|200gT ) ran additional models with 
metallicity enhanced 5 times relative to solar (at odds with 
the measured sub-stellar metallicity (|Torres et al.ll2008l )') or 
a rotation period twice that of th e orbital period (a t odds 
with tidal synchronization theory (R asio et al.lll99g . e.g.)). 
Our fits do not rely on these mechanisms, but rather an 
additional opacity source most likely associated with con- 
densates. 

In an attempt to understand the discrepancy be- 
tween their model predictions and the data at 4.5/im 
IKnutson etafl (|2012T i suggest that advection of CO from 
the dayside to the nightside (where the Carbon should con- 
vert to C H4 in chemical equilibrium ) is i mportant. Sim- 
ulatio ns of ICooper fc Showman! |2006) and lAgundez et all 
j2012T ). studying the carbon distribution, support this claim. 
Give the agreement between our simulations and the obser- 
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vations at 4.5/im on the nightside, it appears that we do 
not require increased CO on the nightside. Moreover, the 
decreased advective efficiencies we find in our simulations 
would self-consistently predict less CO on the nightside as 
well. At the important shorter wavelengths of 3.6/xm and 
4.5/^m, we o ver-predict day s ide fl uxes in comparison to ob- 
servations of lKnutson et al.l (|2012h . If the amount of CH± on 
the dayside was increased above chemical equilibrium values, 
as pr edicted by photochemical models l|Visscher fc Mosea 
l201ll ). the flux at 3.6/^m would be suppress ed bringing our 
model s into better agreement. As noted in iKnutson et all 
l|2012h this would also decrease the flux at 8 and 24^im. 
However, it is important to remember that the addition of 
condensates may mask many details associated with non- 
equilibrium molecular distributions. 

From our analysis of HD 189733b, it is clear that 
the hazes and/or grains play a significant, pe rhaps dom- 
inant, role in determining the overall opacity. IPont et aD 
(2012) suggest molecular opacities may not be relevant for 
HD189733b, and almost all the features of their data can 
be explained simply by a combination of a cloud deck and 
Rayleigh scattering dust grains. Thus it is prudent to con- 
sider if the atmosphere can form grains and hazes. Ideally 
one would self-consistently calculate the formation, destruc- 
tion, advection, and settling of grains and hazes. This is not 
currently included in our model, but is an important direc- 
tion for future research. The vertical velocities in the upper 
radiative zones are relatively modest compared to the zonal 
flows, but when coupled to the large scale-height of the at- 
mosphere, could lead to significant vertical mixing required 
to keep grains at low pressures. 
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